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Introduction 


The gamma-ray spectrum simulation program BSIMUL was designed to allow the operator to 
follow the path of a gamma-ray through a detector, shield and collimator whose dimensions are 
entered by the operator. It can also be used to simulate spectra that would be generated by a 
detector. Several improvements have been made to the program within the last few months. The 
detector, shield and collimator dimensions can now be entered through an interactive menu whose 
options are discussed below. In addition, spectra containing more than one gamma-ray energy can 
now be generated with the menu - for isotopes listed in the program. Adding isotopes to the main 
routine is also quite easy. Subroutines have been added to enable the operator to specify the 
material and dimensions of a collimator. 

This report details the progress made in simulating gamma-ray spectra for a variety of user- 
specified detector designs./ In addition, a short discussion of work done in the related areas of 
pulse shape analysis and the spectral analysis is included. The pulse shape analysis and spectral 
analysis work is being performed pursuant to the requirements of contract F-94-C-0006, for the 
Advanced Research Projects Agency and the U.S. Air Force. 
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The Simulation Program 


A random number generator is used to simulate the behavior of gamma-rays as they pass 
through a detector. The absorption coefficients of the materials making up the detector, shield and 
collimator are used to calculate the reaction probability for: 1) photoabsorption, 2) Compton 
absorption, 3) pair production and 4) Bremstrahlung. The direction that the gamma-ray takes after 
a Compton interaction is calculated using the random number generator in conjunction with 
Equation 1 [Ref 1], while the energy of the gamma-ray is calculated using the Klein-Nishina 


cos0 = 1 +511 keV 


1 

~E 


1 


formula. BSIMUL was originally written with the assumption that the detector and shield would 
be cylindrically symmetric. Taking advantage of the symmetry increased the speed greatly. The 
program can now generate a spectrum resulting from one million gamma-rays entering the detector 
or shield on a SUN Sparc20 in 25 minutes, so speed is no longer critical. New subroutines have 
been added to model a collimator with a hole pattern specified by the operator. While the size and 
locations of the collimator holes must now be entered into a subroutine by the operator, in the 
future they will be entered through the menu. The collimator subroutines may be expanded to 
permit modeling of completely asymmetric detectors. 

Materials with 
absorption constants 
written into the program 
are aluminum, tungsten, 
germanium, silicon, lead, 
sodium iodide, cesium 
iodide, BGO and 
beryllium. The 
absorption constants for 
each material are entered 
for an energy range of 20 
keV to 10 MeV. 

Expanding the program to 
allow the operator to 
simulate detectors and 
shields made of any other 
material requires only a 
knowledge of the 
absorption constants of 
the material. 


Window 



Collimator 


Figure 1 General detector and shield configuration. 
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SIMUL Menu 


The first screen is a table of contents that allows the operator to enter the general information, 
graphics information, and detector and shield configuration screens listed below: 

1 General information for the calculation 

2 Graphics information 

3 Detector geometry 

4 Case and shield geometry 

5 Cold finger geometry 

The second screen, shown in Table 1, contains the number of events to be simulated, the 
isotopes and their relative strengths, and the source location. A typical spectrum contains between 
one million and five million gamma-rays, only a few of which deposit energy in the detector. A 
spectrum of one hundred thousand gamma-rays is useful only for studying the photopeaks. The energy 
entered in line three is used only if a single energy spectrum is desired. If isotopes are specified, the 
energy on line 3 is ignored. Up to ten isotopes may be entered in the menu. The x and y values are 
perpendicular to the axis of the detector shown in Figure 1. z is along the axis. 


Table 1. General information for the simulation 


Number of events 

1000000 

Random seeds 

5676654, 9612341 

Energy (keV) 

1000. 

Number of isotopes 

2 

List isotope and relative intensity 

Co60, 0.5 


Csl37, 0.5 

Initial x and y values 

0.0, 0.0 

Initial z value 

100.0 

Lower limit for BGO detection 

50.00 

Lowest energy for bremstrahlung 

5000.0 

Effective charge for bremstrahlung 

70.0 
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Graphics information, shown in Table 2, is entered on the next screen. The graphics 
option on line 1 allowing the operator to have the path of each gamma-ray plotted as it 
passes through the shield and detector. In this case, the position, momentum, and energy 
of the gamma-ray are printed at each new reaction site. The type of reaction that occurs 
is also printed to standard output. Two histogram files are always generated, the first 
containing all counts and the second containing only counts that were not vetoed by the 
shield or collimator. Gamma-ray energies are converted to channel numbers with the 
calibration information on lines three and four before the histograms are written. 


Table 2. Graphics information 


Graph the events in the detector? 

N 

Histogram file names 

testl.sp, testis. sp 

Energy range for histogram 

20., 3020. 

Size of histogram 

4096 


Table 3 shows the detector material and geometry. At the moment, the detector must be 
cylindrical and made of Ge, Si, Nal or Csl, although other materials can be added easily. 

The inner radius is always zero while the position of the bottom of the detector material is in 
relation to the bottom of the shield, which is always at z=0. The dead layers are only on the 
side and top of the material and may be set equal to zero. All of the lengths are in 
centimeters. 


Table 3. Detector geometry 


Detector material (Ge, Si) 

Ge 

Inner, outer radii of material 

0.0, 2.95 

Bottom, top of material 

4.70, 10.14 

Thickness of dead layer on side 

0.00003 

Thickness of dead layer on top 

0.00003 


Table 4 contains the case, shield and collimator information. Any of these can be 
omitted, but the bottom of the shield must always be at zero with the case above it. Also, 
the inner radius of the shield must be larger than the outer radius of the case. The shield 
can extend any distance beyond the detector. While the collimator material and dimensions 
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are entered here, the hole pattern must be entered in the program itself. 

Table 4. Aluminum case and shield configurations. Collimator information. 


Is there an A1 case? (Y/N) 

Y 

Inner, outer radii of A1 

3.35, 3.51 

Bottom, top of A1 base 

2.01, 2.51 

Length of case 

6.85 

Is there a shield (Y/N) 

Y 

Shield material (BGO, Plastic, etc.) J 

BGO 

Inner, outer radii of shield 

4.00, 6.00 

Bottom, top of shield base 

0.0, 2.0 

Length of shield 

14.0 

Collimator material (BGO, Csl) 

NONE 

Inner, outer radii of collimator 

0.0, 3.515 

Bottom, top of collimator 

21.0, 31.0 


The last screen, shown in Table 5, contains information on the cold finger and window. 

A thin beryllium window has been found to have very little effect on the spectra. The inner 
radius of the cold finger is always zero. 


Table 5. Cold finger and window 































Calculations 


The program will now generate a spectrum with gamma-ray energies from up to ten 
isotopes, the identity and relative strengths of the isotopes being entered by the operator. 

For each isotope, the relative probability of a gamma-ray occurring at a particular energy is 
entered in the program. The energy of each gamma-ray in the simulation is chosen with a 
weighted random number. Isotopes that are not included in the program may be added 
easily in the main routine orsimul.f. The original method of running the program - with one 
gamma-ray energy - can still be used to generate spectra for peak-to-Compton ratio 
calculations. With a 30% HPGe detector as a model, several calculations with C06O and 
Csl37 and BGO and Csl shielding were made. Figure A1 shows a comparison of real and 
simulated spectra with Cs 137 and Co 60 sources for the 30% HPGe detector with no 
shielding. The figures referred to in this section may be found in Appendix A. The 
simulated spectrum was normalized so that the 1 173 keV peak has the same number of 
counts in both spectra. Co60 and Csl 37 calibration sources were used for the real spectrum, 
which was taken with a 30% HPGe detector. The dimensions of the detector, shown in 
Table 1, were used in the simulations. The agreement between the two spectra is excellent 
in the region of 700 keV to 1200 keV. At low energies, the natural background and 
difference in number of Cs 137 gamma-rays causes the spectra differ. Small peaks seen in 
the simulated spectra are the 51 1 keV peak and the first and second escape peaks for the two 
gamma-rays of Co 60. The peaks in the simulated spectrum are at the moment one channel 
wide, but a cubic spline fit to a peak taken with the detector being simulated can be used to 
generate a more realistic peak shape. 

Figures A2 through A5 show the effect of the shield material and thickness on Compton 
suppression. The configuration of the germanium detector used in the calculations is shown 
in Table 7. Figures A2 through A5 were generated by dividing the non-Compton suppressed 
counts by the total counts in each channel of the SIMUL output histograms. A Csl shield 
was used for the calculations shown in Figures A2 and A3. The shield was 1.5 cm thick in 
calculations shown in Figure A2 and 2.5 cm thick in the calculations shown in Figure A3. 
The shield material used to generate the next two figures was BGO, with 1.5 cm thick shield 
used in Figure A4 and a 2 cm thick shield used in Figure A5. 

Also interesting is the ability of a shield to veto gamma-rays coming from the wrong 
direction. FigureA6 shows non-vetoed spectra arising from sources at two positions. The 
detector is the 30% germanium detector used above, with a cesium iodide shield 2.48 cm 
thick and extending 9.86 cm in front of the detector. The larger spectrum was generated 
with the gamma-rays at an initial position of x = 10 cm, y = 0 cm and z = 100 cm (the 
origin is at the base of the shield, so the top of the detector is at 10.14 cm and the top of the 
shield is at 20.00 cm). Gamma-rays entered at an angle of 5.7 degrees from the axis 
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Table 6: Dimensions of the 30% HPGe detector used in the experiments and simulations. 



Minimum value (cm) 

Maximum value (cm) 

Detector: inner and 
outer radii 

0 

2.95 

Detector: bottom and 
top 

4.7 

10.14 

Dead layer thickness 


0.00003 

A1 case: inner and 
outer radii 

3.35 

3.51 

A1 case: bottom and 
top of the base 

2.01 

2.51 

A1 case: length 


7.85 

Shield: Bottom and 
top of the base 

0 

2 

Shield: top 


14 

Cold finger: Inner and 
outer radii 

0 

0.5 

Cold finger: Bottom 
and top 

0 

6.7 


and reached the detector without entering the shield. In the second spectrum, the initial 
positions were x = 50 cm, y = 0 cm and z = 100 cm, so the gamma-rays entered the detector 
at an angle of 26.6 degrees from the axis and had to pass through the shield, which vetoed 
80% of the total counts. 

Figures A7 through A9 show the effect of the length of the shield on the spectra. Figures 
A7 and A8 show the total and non-vetoed spectra for the 30% HPGe detector and BGO 
shield configuration listed in the menu shown on pages two and three. In Figure A7, the 
shield extended four centimeters in front of the detector, while in Figure A8 it extended ten 
centimeters in front of the detector. The total spectra are nearly identical in the two figures, 
but the Compton edges of the non-vetoed spectra are much narrower in Figure A8 because 
the longer shield vetoes gamma-rays that are scattered. A gamma-ray entering the center of 
the detector can be vetoed by the shorter shield if it is scattered at an angle of less than 156 
degrees, while the longer shield can veto a gamma-ray scattered at an angle of up to 170 
degrees. The Co 60 Compton edges of the two non-vetoed spectra are compared in Figure 
A9. 
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Plans 


There remain four major improvements to be made to the program: 

1) Include an option to allow the operator to use cubic splines to generate simulated 
histograms that have accurate peak shapes. At the moment the photopeak lines are one 
channel wide, which is adequate for the calculation of peak-to-Compton ratios and 
germanium spectra, but would not provide a useful simulation of a Csl scintillation 
spectrum, for instance. Several sample spectra have been taken with germanium, sodium 
iodide, and cesium iodide detectors that can be used to test the peak shape subroutines. 
Estimated time: 2 weeks 

2) There are many places where the program could be greatly simplified. In addition, the 
speed may be increased significantly. This is a relatively low priority since it runs in 25 
minutes on a SUN Sparc20 with an input of 1 million gamma-rays and a cylindrically 
symmetric configuration. Estimated time: 4 weeks 

3) The range of materials whose behavior can be simulated can be expanded to include 
CdTe and Hgl 2 , for instance. The photo-absorption coefficients for these materials are 
available and the Compton (a) and pair-production (k) coefficients may be approximated 
from known coefficients as follows: 
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where p is the density, A is the atomic weight and Z is the atomic number. The subscripts 
1 and 2 refer to any two elements. For materials whose photoabsorption coefficients are not 
available, the interpolation formula 
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may be used. [Ref. 2] A study of the accuracy of the interpolation formulas will be 
performed, and the accuracy will be compared to the accuracy of the values of the 
absorption coefficients in use in the program. In addition, a commercially available program 
that calculates absorption coefficients may be used in conjunction with BSIMUL to allow the 
operator to simulate the behavior or any detector or shield material. 

4) Finally, a graphical interface that allows the operator to draw the detector would be a 
great improvement over the current menu system. Popup menus would allow the operator to 
choose the material of each part and the exact dimensions would be entered through a text 
window. 

In order to make further progress, we need the detector configurations for which spectra 
are desired. In addition, we need to know which materials are being consideied for detectors 
and shielding. A list or sources and their probable strengths, or the environments that the 
detectors will be used in will enable us to perform realistic simulations. 
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Progress on Related Work 


Pulse shape analysis 


In a typical gamma-ray spectrum, one finds several undesirable features. In addition 
to the desired photopeaks, there are a Compton continuum background and Compton edge 
features. These features can be considered noise. As in any information bearing medium, 
one would benefit from increasing the signal to noise ratio. In a gamma-ray spectrum, this 
could be accomplished by reducing the Compton continuum and the Compton edges. 

A newly developed method of doing this is pulse shape analysis. A pulse from a 
gamma-ray detector may contribute to the photopeak, the Compton continuum or the 
Compton edges. Rejecting pulses which contribute to the Compton continuum or Compton 
edges would eliminate the noise of the spectrum. By applying pulse shape analysis to each 
pulse coming from the detector, one may reject some of the undesired pulses and effectively 
improve the signal to noise ratio of a given spectrum. 

Pulse shape analysis can be used to reduce the peak to noise ratio in a spectrum by 
vetoing those pulses that come from the edges of the detector. Gamma-rays that interact 
with the detector near the edge are much more likely to deposit only part of their energy in 
the detector than those that interact with the detector near the center. We have performed a 
series of experiments to optimize the criteria for determining the location of a pulse in a 
detector. 

One criterion is the rise time. Pulses caused by interactions which occur near the 
outer edge of a detector (N-type cylindrical HPGe) tend to have a longer rise time than 
interactions occurring near the center of the detector. Therefore, rejecting pulses with long 
rise times will reject photons which have the greatest chance of scattering out of the 
detector. The method has been demonstrated to work quite well, with Peak to Compton 
ratios improving by a factor of 1.5. Other features of the pulses that we examined were the 
fall time and the widths of the pulse at 70% of the maximum and at 50% of the maximum. 
These criteria are almost completely energy independant. The ratio of the 70% pulse 
maximum to the 50% pulse maximum increases at the photopeaks and decreases at the 
Compton edges. The Compton features of a spectrum can be reduced by using these criteria 
to veto pulses. Figure 2 shows an unsuppressed spectrum and an unsealed suppressed 
spectrum. Figure 3 shows the same two spectra, only the suppressed spectrum has been 
scaled so that the 662 keV photopeaks of the two spectra are equal in height. 

Presently at Constellation Technology, other, more complex criteria are being 
developed which promise even higher peak-to-Compton ratios. Some of these have been 
demonstrated in digital experiments and have been implemented in analog hardware, 
allowing real time analysis. 
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Figure 2. Upper Curve: Upsuppressed Spectrum 

Lower Curve: Suppressed Spectrum 

Sources: C 06 O and Csl37 
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Figure 3. Upper Curve: Unsuppressed Spectrum 
Lower Curve: Suppressed Spectrum scaled so that 
the Csl37 Peaks are the Same Height. 




ROBFIT 


The robust spectral fitting program ROBFIT uses cubic splines to fit gamma-ray 
spectra down to the noise level. The program is separated into two main routines: the 
background and peak fitting routines. Cubic splines are used to fit the background, while up 
to ten standard peak shapes are used to fit the peaks. The available peak shapes are 
Gaussian, Lorentzian, Voight and a computer-calculated peak shape that uses cubic splines to 
fit a peak in the spectrum that is specified by the user. The fits are performed as follows. 

1) Calculate the standard peak shapes. 

2) Fit the background with up to sixty constants that specify the cubic spline positions and 
magnitudes. The background function can be written [Ref 3]. 

3 m 4 

S(x) = £ ex' + Y,d.(.k.-xf 

1=0 ' y=l 1 } 


where the number of splines, M, is specified by the user. 

3) Subtract the background fit from the spectrum and fit the peaks with the standard peak 
shapes. 

4) Subtract the fitted peaks from the spectrum and refit the background. 

5) Repeat the last two steps until the error is within the limit set by the user. 

ROBFIT has been used successfully in the analysis of gamma-ray data from the 
supernova SN 1987 A [Ref 4], The spectrum shown is Figures 4 and 5 show a typical 
spectrum and fit from the supernova [Ref 3]. 

Work on the spectral analysis program ROBFIT has proceeded along two lines. 

First, a new fitting program that compresses the data has increased the speed of the program 
by a factor of 10. Second, a new user interface has been written that makes the program 
much easier to use. Modifications to the data display routines now allow the operator to 
change scales with the arrow and page keys on both UNIX and DOS machines. It is now 
possible to use cubic splines to calculate the full width-half max of a peak while running the 
raw data display program. 

The data compression routines are used by both the background and peak fitting 
portions of the program to increase both the speed and the accuracy of the fits. The original 
program compressed the background linearly to 1024 points, while the new routines perform 
a non-linear compression that places more points in areas where the background varies 
quickly and fewer points where the background is relatively constant. A new option added 
to the background fitting program is FBKG, which fits the background of a spectrum that is 
flat at large channel number. The spectrum shown in Figure 6 has been fitted with forty 
background constants and a cutoff of 100 by both versions of the background program. The 
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Figure 4. Gamma-Ray Spectrum From Supernova SN 1987A [REF 3] 
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Figure 5. Fit Over a Small Range of the Supernova Spectrum [REF 3] 


program ran in 13 minutes with non-linear data compression and in 26 minutes with linear 
compression. It has been found that compressing the peaks also results in an increase in 
speed without loss of accuracy. The compression factors for both the peak and background 
programs will be optimized. 
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Background Spectrum 



Figure 6. Background fit to the compressed data, 
constants were used. The cutoff was 100. 
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Appendix A 


Simulated Co 60 and Cs 137 Spectra 
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Real and Simulated Spectra 



Figure A1. Curve 1: simulated spectrum normalized at the 
1173 peak. Curve 2: real spectrum taken with 30% HPGe 
detector. Sources: Co60, Cs137 














Noise Reduction Due to Shielding 






re A6. The upper cu 
axis. The lower cur' 
axis. Sources: C06O 












